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І. INTRODUCTION 


A. BACKGROUND 


The ability to obtain localization information from a cellular phone transmission 
is increasingly becoming a much sought-after facility. Thanks to some improved 
techniques in cellular communications, this ability is very much realizable. As a matter 
of fact, recently there has been legislation which requires the development of this 
technology. When a 911 call is made from a wire line phone, the dispatcher 
automatically can provide the location of where the call was made. Currently, there is no 
such capability for calls made from cellular phones. In October 1999, President Clinton 
signed a law which makes 911 the official nationwide emergency number for both 
regular and cellular phones [Appendix A]. This legislation followed a move by the FCC 
which requires, by the end of 2001, that cellular 911 calls automatically provide the 
caller's location [Appendix A]. 

There are countless ways in which the development of this technology will be 
beneficial to many different public and private enterprises. If this technology is utilized 
to 1ts full potential, 1t will completely change the way the entire wireless industry does 
business as well provide an invaluable resource for many military applications. Listed 
here are just a few of the benefits: 

а Location Sensitive Billing — Network operators will be able apply charges 

depending on the location of the mobile. In this way, different rates can 


be charged for calls from the home or office. This would allow network 


operators without copper-based public switched telephone network 
(PSTN) to offer competitive rates for calls from the home or office.[1] 
Intelligent Transport Systems — Traffic congestion can be determined by 
the number of cellular calls coming from a particular stretch of highway. 
Vehicles with installed maps can provide directions and alternate routes to 
travelers to avoid congestion and traffic mishaps. Also companies could 
track the location of all the vehicles in its fleet and more efficiently 
dispatch vehicles according to the traffic situation and distance to the final 
destination. 

Enhanced Network Performance — By monitoring the frequency of calls 
from various locations over time, cellular networks can be better 
organized to provide the optimum service for a particular coverage area. 
This includes the ability to accurately monitor the movement of the mobile 
telephone to enable the network to make better decisions on when to 
initiate call handoff to a different cell. 

Electronic Surveillance — Law enforcement authorities can use position 
information to track the location of a mobile station which could serve a 
multitude of purposes.[2] 

Automobile Security — Location information can be gained from a 
strategically placed mobile transmitter in the event of a car theft. This 
would aid tremendously to the recovery of stolen vehicles and greater rate 


of apprehension of car thieves. 


в Unauthorized Cellular Usage Detection – Cellular carriers have long been 
plagued by the unauthorized access to personal communications systems . 
These illegal activities would diminish drastically with the advent of 


cellular phone localization. 


B. OBJECTIVE 


There are several methods available that would allow us to implement a position 
location capability into cellular communications systems. For reasons that will become 
more apparent in the following chapters, the Time Difference of Arrival (TDOA) method 
appears to be the most logical choice. In this thesis, we will discuss the TDOA 
estimation and its application to Global System for Mobile (GSM) signals. Specifically, 
we will investigate several methods to reduce the mean squared error of the TDOA 
estimate, thereby allowing a position localization with a corresponding lower error 
probablity. 

GSM was chosen to provide the backbone for this research because of its 
worldwide acceptance as the Personal Communication Service (PCS) network of choice 
and its inherent positioning localization advantages. We will discuss this system in great 
detail in order to provide an understanding of its advantages. We will also demonstrate 
our ability to simulate GSM signals using the Hewlett-Packard Advanced Design System 
(HP-ADS) program. Finally, these generated GSM signals will be evaluated using 
Matlab-based programs to provide an indication of how to improve upon the localization 


of a GSM emitter. 


E: RELATED WORK 


The initial research on this project was begun in 1998 by Unal Aktas, LTJG, 
Turkish Navy. His work focused on the use of the wavelet transform to increase the 
accuracy of the TDOA estimation. He investigated several different denoising techniques 
and presented a comparison of the techniques. These included wavelet denoising based 
on рана method, wavelet denoising using the Wo-So-Ching threshold, wavelet 
denoising using hyperbolic shrinkage, wavelet denoising using median filtering, 
modified approximate maximum-likelihood delay estimation, denoising based on the 
fourth order moment, and a time-varying approximate maximum-likelihood technique 
The results of his work showed that the fourth order moment, and the modified 
approximate maximum-likelihood techniques were superior to the others. For that 
reason, this thesis only makes use of those three techniques. Modified versions of the 


Matlab code for wavelet-based denoising developed by Aktas, are used in this thesis. 


П. GSM AND THE TIME DIFFERENCE OF ARRIVAL 


A. REVIEW OF GSM 


Global System for Mobile Communications (GSM) is the European standard for 
digital mobile telephones. It was developed in 1990 and first implemented in 1991 to 
replace first generation European cellular systems, which were generally incompatible 
with each other. It has gained worldwide acceptance as the first universal digital cellular 
system with modern network features extended to each mobile user, and is a strong 
contender for Personal Communication Service (PCS) above 1800 MHz throughout the 


world. [3] 


1. GSM Channel Structure 

One characteristic of GSM that makes it so desirable is that it uses both Time 
Division Multiple Access (TDMA) and Frequency Division Multiple Access (FDMA) to 
transmit and receive information. TDMA allows for multiple conversations to take place 
simultaneously at the same frequency (channel) by rationing different time slots to each 
user. GSM uses two 25 MHz bands of the radio frequency spectrum, one for uplink 
(mobile station to base station) and the other for downlink (base station to mobile). 
FDMA divides each of these blocks into 125 channels, each with a bandwidth of 200 
kHz. TDMA allows for 8 different time slots , as shown in Figure 2.1 from Ref [1], 
with each time slot being 576.92 microseconds in duration. The grouping of these 8 time 


slots 1s called a frame. This combined with the mechanics of FDMA provides for a total 


of 1000 actual speech or data channels. Currently, there is a half-rate speech coder in the 


process of development that will double the number of channels to 2000 [4]. 
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Figure 2.1: GSM Frame Structure 

Of the total number of available channels, a certain number of them are dedicated 
as control channels or logical channels. Each logical channel has a specific purpose 
ranging from transmitting system parameters, to coordinating channel usage between 
base station and mobile. The logical channels are mapped onto dedicated time slots of 
particular frames. This is executed using a burst structure, the most commonly used 
being a normal burst. Contained in the middle of this burst is a 26-bit training sequence, 
chosen for its correlation properties. This received burst is correlated with a local copy of 
the training sequence to allow for the estimation of the impulse response of the radio 


channel, which aids in the demodulation of the bits in the burst. Another important 


characteristic of training sequences in the different burst structures of GSM is that they 


lend themselves to the measurement of time which can be used to aid in the gathering of 


time-based positioning information [1]. 


2 GMSK 


GSM uses 0.3 Gaussian Minimum Shift Keying (GMSK) as its method of 
modulation. This facilitates the use of narrow bandwidth and coherent detection 
capability. In GMSK, the rectangular pulses are passed through a Gaussian filter prior to 
their passing through a modulator. Within this filter, baseband Gaussian pulse shaping 
smoothes the phase trajectory of the MSK signal, and thus stabilizes the instantaneous 
frequency variations over time. The result of this is significantly reduced sidelobe levels 
in the transmitted spectrum. The value 0.3 represents the 3 dB-bandwidth-bit duration 
product (BT) or normalized pre-Gaussian bandwidth, which corresponds to a filter 
bandwidth of 81.25 kHz for an aggregate data rate of 270.8 Kbps. As shown in Figure 
2.2 from Ref. [4], sidelobe levels fall off dramatically with decreasing BT product. 
However, reducing the BT product increases the irreducible error rate which is a product 
of intersymbol interference (ISI) created in the Gaussian filtering process. GMSK 
sacrifices the irreducible error rate for extremely good spectral efficiency and constant 
envelope properties. This is a good tradeoff as long as the irreducible error rate is less 
than that produced by the mobile channel. 

All base stations transmit a GSM synchronization burst on the synchronization 
channel. Contained in the burst is another useful 64-bit training sequence. The 
autocorrelation function for this sequence reveals a main correlation peak with a width of 


4 bits. Compared to signals used in other time-based positioning systems, this is 


relatively wide. The GMSK modulation technique 1s responsible for this increased width. 
This is significant from a time-based positioning standpoint because it reduces the timing 


resolution and hence the maximum achievable position accuracy. 
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Figure 2.2: Power Spectral Density of aGMSK Signal 


3. Interference 


In the mobile radio path, the GSM signal will experience some ISI. The mobile 
radio channel is generally hostile in nature and can affect the signal to such an extent that 
the system performance is seriously degraded [5]. Inherent in any communication system 
is the phenomenon known as multipath. This refers to the situation where a signal 
propagating from a transmitter to a receiver can travel via several different paths 
including line of sight and one or more reflected paths. The reflected paths cause delay, 
attenuation in the amplitude, and phase shifts relative to the direct path. In some 
environments, multipath components may be delayed by up to 30 microseconds. These 
delayed components cause ISI in the received signal. [6] An example of this interference 
1s illustrated in Figure 2.3 from Ref. [1]. As shown, it is a distorted correlation peak 
which can cause errors in any time measurements. The degree of distortion depends on 
the relative amplitude, delay, and phase of the received signals. 

The conventional GSM receiver contains a signal demodulator which equalizes 
the received signal to improve system performance in the presence of multipath by 
combining the information received from the different arrivals. This process 1s called 
adaptive equalization. The equalizer does not reject multipath signals, but combines 
them. However, for positioning reasons, the receiver needs to be able to determine the 
first arrival corresponding to the line of sight, and thus requires a multipath rejection 


algorithm [1 ]. 
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Synchronization Burst (64-Bit Training Sequence) 


4, GSM Coverage Area Layout 


The area that a GSM network covers is divided into a number of cells, each 
served by its own base station (BS), also known as a base transceiver station. The 
number of cells served by each BS, called a cluster, is determined by the amount of 
system capacity (total available channels) required for a particular area. A larger cluster 
size is an indication of a larger capacity system. One tradeoff associated with the larger 
capacity system is that a more complex system is required to control the increased 


demand for frequency reuse. Also, more cells per cluster translates into smaller 


individual cell size and thus, more co-channel interference for which to account, as well 
as a greater frequency of handoffs between cells. The typical cluster size 1s 4, 7, or 12. 
Each cell is distinguished by a unique cell identifier and is allocated one or more 
uplink/downlink frequency pairs. Multiple BS’s may be grouped together under the 
control of a base station controller (BSC). In turn, several BSC's are usually controlled 
by a mobile service center (MSC), which handles tasks such as call routing and serves as 


the interface between the mobile network and the fixed telephone network. 


Э: GSM Performance Enhancing Techniques 

In addition to increasing the number of cells per cluster, GSM uses a number of 
other methods to increase capacity. These are important because they provide favorable 
implications for positioning. One such technique is the use of sectored cells. Sectoring 
effectively breaks the cell down into groups of frequencies (channels) that are only used 
within a particular region. This is best exemplified as slicing a pie into equal segments. 
Not only does this technique offer additional information for use by a positioning system, 
but it also has the effect of reducing co-channel interference. This reduced co-channel 
interference is due to the fact that a given cell receives interference and transmits with 
only a fraction of the available co-channel cells. 

Another performance enhancing technique employed by GSM is slow frequency 
hopping; so termed because the hopping rate 1s low compared with the symbol or bit rate. 
The mobile radio channel is a frequency-selective fading channel, which means the 
propagation conditions are different for each individual radio frequency. Slow-frequency 


hopping is used instead of fast-frequency hopping because, in order to perform well, a 
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frequency synthesizer must be able to change its frequency and settle quietly on a new 
one within a fraction of one time slot (576.92 microseconds). Frequency hopping 
provides frequency diversity to overcome Rayleigh fading due to the multipath 
propagation environment. Frequency hopping also provides interference diversity. The 
amount of interference varies with any given channel within any given cell. A receiver 
set to a channel with strong interference will suffer excessive errors over long strings of 
bursts. Frequency hopping prevents the receiver from spending successive bursts on the 
same high-interference channel. These beneficial effects of frequency hopping 
contribute to a reduction in the minimum signal-to-noise ratio (SNR) required for good 


communications. 


B. POSITIONING TECHNIQUES 


It is important to remember that in order to comply with the FCC ruling, any 
positioning solution must work with existing phones. For this reason, most cellular 
geolocation proposals are multilateral, where the estimate of the mobile’s position is 
formed by a network rather than by the mobile itself. Hence, the unilateral proposal of 
adding a global positioning system (GPS) capability to mobile phones is not a universal 
solution since mobile telephone operators would have to replace or retrofit every mobile 
telephone. Even if it were feasible to carry out this retrofit, there are other factors that 
would make this a futile attempt. First of all, the size and cost of mobile transmitters 
would be considerably increased because of the incorporated GPS receiver. Moreover, a 


clear view of the sky is necessary to receive usable GPS data from at least three satellites, 
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which makes the system useless in areas surrounded by buildings, mountains, and other 
obstructions—common working environments for cellular communications. [7] 

There are three ways that position localization information can be determined 
using the signaling aspects inherent in the GSM specification. One is through the Angle 
of Arrival (AOA) method, which uses sector information. Second, is through the use of 
propagation time using timing advances (TA). Third, is the Time Difference of Arrival 
(TDOA) method, which offers several advantages over the previous methods. In addition 
to providing some immunity against timing errors when the source of major reflections is 
near the mobile, the TDOA method is less expensive to put in place than the AOA 
method. Also, the AOA and TOA methods are hampered by requirements for line of 
sight (LOS) signal components, whereas the TDOA method may work accurately without 
a LOS component [8]. Each method will be briefly described before taking an in-depth 


look at the TDOA 


1. Angle of Arrival (AOA) 


If two or more BS platforms are used, the location of the desired mobile in two 
dimensions can be determined from the intersection of two or more lines of bearing 
(LOB), with each LOB being formed by a radial from a base station to the mobile. 
Figure 2.4 illustrates the method known as triangulation using three different platforms. 
Here, the LOB”s from the mobile to two adjacent platforms form what is called a 
measured angle, y. After two measured angles are calculated, trigonometry or analytic 


geometry may then be used to deduce the location of the MS. 
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Figure 2.4: Triangulation 

Another variation of the AOA method relies upon the sectoring of cells for all of 
its positioning information. Cell sectoring involves the use of highly directive antennas 
to effectively divide the cell into multiple sections. The angle of arrival can be 
determined at the base station by electronically steering the main lobe of an adaptive 
phased array antenna in the direciffh of the arriving mobile signal. In this case, a single 
platform may be sufficient for AOA positioning, although typically, two closely spaced 
antenna arrays are used to dither about the exact direction of the peak incoming energy to 
provide a higher resolution measurement. [9] 

AOA measurements might provide an excellent solution for wireless transmitter 
localization were it not for the fact that this method requires signals coming from the 


mobile to the base station be from the line of sight (LOS) direction only. This is seldom 
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the case ın cellular systems given that the operating environment is often heavily 


cluttered with rough terrain, tall buildings and other obstructions. 


2 Time of Arrival (TOA) 

Since GSM is a TDMA system, its successful operation depends on the ability of 
all signals to arrive at the BS at the appropriate time. And since the signals arriving at the 
BS’s originate from different distances, the time at which the signals are sent must be 
varied. This is accomplished by having each BS send a timing advance (TA) to each MS 
connected to ıt. The TA is the amount by which the mobile station (MS) must advance 
the timing of its transmission to ensure that it arrives in the correct time slot. Obviously, 
the TA is a measure of the propagation time, which 1s proportional to the distance from 
the BS to the MS. Hence, the location of the MS can be constrained to a circular locus 
centered on the BS. If the MS could somehow be forced to hand over to two more BS's, 
it would be possible to implement a positioning scheme based on the intersection of the 
three circles. 

One problem with the TOA technique is that artificially forced handoffs can be 
made to less optimal BS's. This degrades call quality as well as reduces system capacity. 
Also, this method requires that all transmitters and receivers in the system have precisely 
synchronized clocks, since just 1 microsecond of timing error could result in a 300 meter 
position location error [9]. The use of the TA, which is in essence a timestamp, is a 
burden that most would rather not have to deal with. Furthermore, the employment of it 
presents another problem. Under GSM specifications, the TA 15 reported in units of a bit 


period, which equates to a locus accuracy of 554 meters. Even this considerable amount 
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of area will increase when multipath degradations and dilution of precision (the amount 


by which errors are degraded by geometry) are taken into account. 


3. Time Difference of Arrival 


The TDOA is basically an improvement upon the TOA method. The idea behind 
TDOA is to determine the relative position of the mobile transmitter by examining the 
difference in time at which the signal arrives at multiple base station receivers, rather 
than the absolute arrival time. Therefore, each TDOA measurement determines that the 
transmitter must lie on a hyperboloid with a constant range difference between the two 


receivers. The equation for this range difference is given by 


R. = A, AM (а, - 22 - (€; -xy +(7, – у) «(Z; -zy. C-D 


where the coordinates (X., Y., Z.) and (ХҮ, 2) represent the position-fixed receivers i 
and j, and x, y, and z represent the unknown coordinate of the target transmitter. [9] 

The location of the mobile can be estimated in two dimensions from the 
intersection of two or more hyperboloids generated using TDOA measurements from 
three or more base stations. This can be extended to three dimensions if four or more 
TDOA measurements are used. 

As in the TOA method, the TDOA requires that all fixed transmitters and 
receivers in the system have precisely synchronized clocks. This, however, corresponds 
to the timing standards already in place at cellular base station sites. Unlike the TOA 


measurements, TDOA measurements do not require a timestamp of the transmitted 


signal. 


16 


IH. POSITION LOCALIZATION USING TDOA ESTIMATION 


The process of localizing the position of a cellular transmission occurs in two 
stages. First, an accurate estimate of the TDOA measured values must be determined 
from noisy signals that are received from the mobile. This involves some type of 
denoising technique, the most popular being wavelet denoising, which is discussed 
thoroughly in [10]. Generally speaking, this involves taking the wavelet transform of the 
received signal, thresholding the wavelet coefficients, and performing the inverse wavelet 
transform on the modified coefficients. Second, a location of position estimate is 
determined from the TDOA calculations. Both phases are very complex and subject to 
minor inconsistencies that can propagate into large, unacceptable errors. Thus, accurate 


and efficient algorithms are required for both stages of processing. 


A. TDOA COMPUTATIONS 


In order to be able to evaluate the hyperbolic range equations, Eq. (2.1), one must 
first estimate the range differences Rij, or equivalently, the TDOA t; — t. The most 
widely accepted method for obtaining these values is the generalized cross-correlation 
method. If the transmitted signal is s(t), then the signal that arrives at one receiver, x(t), 


will be delayed and corrupted by noise and is given by the equation 


x(t) - as(t- di) - nt, 0<о0<1 (3.1) 


ІШ 


where d; is the delay time from the mobile to receiver one, n;(t) is noise, and a is the gain. 


Similarly, the signal that arrives at another receiver, y(t), is given by the equation 
y(t) = Bs(t— dj) + пу), O<B<1 (3.2) 


where d; is the delay time from the mobile to receiver two, n;(t) is the noise, and ß is the 
gain. The cross-correlation function between the two received signals is derived by 
integrating the lag product of the two signals for a sufficiently long time period, T. This 


function is given by 


T 


Ё, (т) = = [x@y@-2)at [9]. (3.3) 


0 


This approach imposes the requirement that the receiving base stations share a precise 
time reference and reference signals, but does not impose any requirement on the signal 
transmitted by the mobile. We can improve the SNR of the TDOA estimation by 
increasing the integration interval, T. The maximum likelihood estimate of the TDOA is 
obtained from the computed cross-correlation function. This estimate is equal to the 


value of t that maximizes the cross-correlation function, Eq. (3.3). 


Another method of estimating the cross correlation function, which yields the 
same results, is to compute the estimated cross-spectral density function. An inverse 


Fourier Transform of the cross-spectral density yields the estimated cross-correlation 


function.[9] 


18 


Prior to the cross-correlation function computation, some frequency domain 
processing may be used to improve the ability to distinguish the desired signal from 
arriving multipath components. This benefit is a result ofthe conduciveness of frequency 
domain processing to filtering techniques which are effective in reducing the effects of 


noise and interference on the TDOA estimate. 


B. OBTAINING A POSITON LOCATION ESTIMATE 


After an accurate estimate of the TDOA measurement is obtained, it is then 
possible to determine the corresponding solution to the hyperbolic equation. This process 
is initiated by substituting the range difference estimate from Eq. (2.1) into the hyperbolic 
equation, Eq. (3.3), then solving for the Cartesian coordinates of the mobile. The 
solution obtained will not be a trivial one, since the range estimates will most likely be 
noisy or inconsistent. The solution to the hyperbolic equation is also complicated by the 
fact that the fixed location receivers are arranged in a non-uniform fashion. This may be 
accounted for by linearizing the equations through Taylor-series expansion in which only 
the first two terms are retained. However, this can result in significant position location 
errors due to geometric dilution of precision (GDOP) effects. GDOP occurs when a 
relatively small range of errors results in a large position location error due to a very large 
difference in distance between the mobile unit location on the hyperbola and the locations 


of the two receivers being used. 
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There have been several proposed solutions to reduce the above inaccuracies. 
One being an exact solution for the case where the number of equations equals the 
number of unknown coordinates. However, this solution does not account for redundant 
information. Other proposals include spherical interpolation, a divide-and-conquer-based 
solution, and more recently, a proposal, by Y. T. Chan, is a closed form solution which is 
valid for both close and distant sources [9]. Chan’s method has been shown to be the 
superior of those mentioned by virtue of a lower computational complexity and higher 
noise threshold tolerance. In fact, Chan’s method guarantees optimum performance 
around small TDOA noise region since 1t has been shown to achieve the Cramer-Rao 
lower bound in this region. [12] 

The essence of Chan's method is to introduce an intermediate variable to the 
original TDOA equations, then transform the nonlinear equations relating TDOA 
estimates and the source position into a set of equations which are linear in the unknown 
parameters and the intermediate variable. A least-squares computation produces an 
initial solution. By exploiting the known relation between the intermediate variable and 
the position coordinates, a second weighted least squares computation gives the final 


solution for the position coordinates. [13] 


e: DETERMINING THE ACCURACY OF POSTION LOCATION 


MEASUREMENT 


In an effort to measure the performance of TDOA position location systems, there 


are several grading criterion at one”s disposal, all of which focus upon the accuracy of the 
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position measurement. The most obvious and most widely used measure of accuracy is 
the root mean square (rms) accuracy. This is defined as the rms deviation of the 
measured position (x,y) about the true position (x,y), and is given by the following 


equation: 


c -4E((*-x) O -»Y] (3.4) 


where E is the expectation. The rms accuracy is a function of the accuracy of the raw 
locus measurements and the geometry of the base stations relative to the mobile. There is 
a limit on the achievable accuracy imposed by errors in the locus measurements. This 
accuracy can be further degraded by the physical geometry of the link between base 
stations to the mobile, referred to earlier as the GDOP. The GDOP for a 2-dimensional 
hyperbolic position location system is defined in [9] by: 


2 2 
+O 


O, y 
ра 7 (3.5) 
OR 


2 2 “+ ° 4 : 5 
where o, and ø, are the mean square position location errors in the x and y directions, 


respectively, and m is the mean square TDOA ranging error. Ideally, it would be 


prudent to minimize GDOP. However, to do this would be infringent upon a very 
important principle in the design of mobile communication systems, where the goal is to 
optimize the availability of service to mobile telephone users. Thus, we can minimize the 
GDOP, but we will degrade the availability of service in the process. Therefore, the 
focus of the accuracy determination must be placed primarily on error in the locus 


measurement. [9] 
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IV. METHODS TO IMPROVE TDOA ESTIMATION 


Any effort to improve the TDOA estimation must focus primarily on reducing the 


effective noise at the receivers. Earlier research presented seven different denoising 


schemes designed toward this end: 1) wavelet denoising based on Donoho’s method, 2) 


wavelet denoising using the Wo-So-Ching threshold, 3) wavelet denoising using 


hyperbolic shrinkage, 4) wavelet denoising using median filtering, 5) modified 


approximate maximum-likelihood delay estimation, 6) denoising based on the fourth 


order moment, and 7) a time-varying approximate maximum-likelihood technique. The 


results of simulation of each of these methods is shown in Figure 4.1. [10] 
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Figure 4.1: Comparison of Denoising Methods for the GSM Signal 
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It was shown that method 1 failed at low SNR’s. Experimental results from [10] 
showed that method 2, while an improvement over method 1, looses its advantage with 
increasing carrier frequency. Additional improvements are made using methods 3 and 4, 
but as can be seen ın Figure 4.1, the best results (lowest mean-square error) are obtained 


using methods 5 through 7. We will take a closer look at these last three methods. 
A. FOURTH-ORDER MOMENT WAVELET DENOISING 


Inevitably, any transmitted signal will acquire some type of additive noise before 
reaching the receiver. However, it is possible to estimate the frequencies of the corrupted 
complex sinusoidal signal. A Gaussian signal being completely characterized by its 
second order statistics and the odd order moments being equal to zero for a symmetric 
probability density function, the separation of the signal and the noise requires at least the 
use of the fourth-order moments [11]. 

To define the fourth-order moment, we first model the received signal as: 

z(n) = х(п) + [(и) (4.1) 
where x(n)ıs a complex zero-mean, non-Gaussian, fourth order stationary signal and 
I(n) is the noise, which is a complex zero mean-Gaussian signal independent of x(n). 


The fourth-order moment is then [11]: 


M, -ük-il-i) s Е(2,,2,,2, ,2,) (4.2) 
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It was shown in earlier research that the fourth order moment of a detail function 
which contains the signal should be greater than 30, * where См ' is the noise power at 


subband d,. [10] 

By using this property, the wavelet coefficients that represent noise can be 
eliminated while those having a signal dependency are retained. After modifying the 
detail functions, the denoised signal is obtained by performing an inverse wavelet 


transform using the modified coefficients. 


B. MODIFIED APPROXIMATE MAXIMUM LIKELIHOOD DELAY 


ESTIMATION 


Critical to the task of source localization is the time delay estimation between 
signals received at two spatially separated sensors in the presence of noise. If we let s(n) 
represent the source signal, nj(n) апа пп) represent the additive noises at the respective 
sensors, and D is the difference in arrival times at the two receivers, then the receiver 
outputs, rj(n) and r2(n), are estimated by: 

г(п) = одз(п) + пи(п), -n=0,1,...,T-1 (4.3) 

r2(n) = asín - D) + nx (mn), 0<04,02<1 (4.4) 
where T is the number of samples collected at each channel. 

In the approximate maximum likelihood delay estimation, after wavelet 
decomposition and prior to cross-correlation, both the channel outputs are optimally 
weighted at different frequency bands with the use of an orthogonal wavelet transform. 


This weighting is done to attenuate the noise spectral components. The scaled subband 


25 


components are then combined using inverse wavelet transform to construct the Modified 
АМІ, (МАМІ) prefiltered signal. The orthogonal wavelet transform is attractive because 
in addition to being computationally efficient, 1t does not suffer from the performance 
degradation due to errors that may occur in the estimation of D used in Eq. (4.4), since in 
practice, a highly accurate estimate of D may be difficult to achieve. [12] 

The final MAML delay estimate 15 calculated from the location of the peak of the 


cross-correlation function of the two denoised signals. 


C: TIME-VARYING MODIFIED APPROXIMATE MAXIMUM 


LIKELIHOOD DELAY ESTIMATION 


The time-varying MAML method, as is obvious from the name, works in much 
the same way as the MAML, except that it additionally accounts for the time-varying 
characteristics of the signal. Due to the time-varying nature of some signals, the 
subbands may not contain actual signal components throughout the entire segment. 
Therefore, the MAML method may loose some of its effectiveness when applied to 
segments in the presence of transitory signals. The time-varying MAML method, in 
essence, allows for the MAML method to be applied twice to each segment, thus 


approximating a time-varying gain. 
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V. SIGNAL SIMULATION AND TEST RESULTS 


A. SIGNAL GENERATION 


All signals used in this project for test purposes were generated using the Hewlett- 
Packard Advanced Design System (HP-ADS). This provided an excellent method of 
obtaining signals that are essentially the same as would be encountered in an actual 


cellular communication receiver. 


1. Advanced Design System 

The HP-ADS program is an invaluable tool for the engineers representing many 
different design aspects such as communications, digital signal processing, electronic 
circuits, mechanical circuits, and many others. Of course, for this project, we were 
particularly interested in the communications package. This package allows the user to 
custom design any type of communications system, run simulations of the design, and 
extract data collected from the simulation. Figure 5.1 shows the system used for the 
signal simulation. A more detailed view of each of the major components of this system 
can be found in Appendix [B]. Appendix [B] also provides detailed instructions for using 


the HP-ADS program to design and simulate a communications system. 


2: GSM Signal Generation 


One parameter that must be specified for the GSM signal generation is the sample 
time. The specifications of the GSM signal were presented in Chapter II. The HP-ADS 


program specifies a symbol period of 3.7037 microseconds. The filter bandwidth is 1.2 


27 


MHz. This allows us to sample at a minimum frequency of 2.4 MHz without violating 


the Nyquist criterion. For test purposes, we use three different sampling frequencies. As 


shown in Table 5.1, these are: 10.8 MHz, which corresponds to a sample time of 92.592 


microseconds; 5.4 MHz, which corresponds to a sample time of 185.185 microseconds; 


and 2.7 MHz, which corresponds to a sample time of 370.370 microseconds. 
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Figure 5.1: HP-ADS GSM Communications System 


28 


Samples/Symbol Sample Time Sampling Frequency Symbols/600 
(microseconds) (Megahertz) Samples 


22592 


185.165 


370.37 





Table 5.1: GSM Signal Parameter Combinations 

The I and Q channel outputs of each of these GSM signals are shown in Figure 
5.2. One fact that is immediately evident from Figure 5.2 is that as the sampling interval 
increases, so does the variation in the GSM signals. The reason for this 1s obvious if we 
refer to Table 5.1. We can see that as the sample interval increases, we can expect to 
contain more symbols within a given number of samples. For test purposes, we use a 
constant number of samples (600) for all simulations. We also realize a better correlation 
characteristic with increasing sample time as shown in Figure 5.3. Better correlation is 
demonstrated by a sharper main lobe with lower sidelobes. Figure 5.3 (d) shows the 
correlation of the Matlab generated GSM signal presented by Aktas [10]. This figure 
represents the best correlation. Plots of the power spectral densities of these signals can 


be found in Appendix [C]. 
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(b) I channel; Ts = 92.592 ns (a) Q channel; Ts = 22.592 ns 





(d) I channel; Ts = 185.185 ns — а Ts = 185.185 ns 
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Figure 5.2: HP-ADS Produced GSM Signals 
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(a) Correlation of sıgnals x and y, ADS 92.59 nsec (b) Corelaton of sgnalsx and y, ADS 185.185 nsec 
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(c) Correlation of signals x and y; ADS 370.37nsec (d) Correlation of signals x and y; MATLAB 
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Figure 5.3: Auto-Correlation Results of Test Signals 


B. SIMULATIONS 


In Chapter IV, we discussed three methods to improve the TDOA estimate. In 
this section, we will use these methods to test our HP-ADS generated signals and see if 


we are able to extract useful TDOA data. Bear in mind that our criterion for 
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improvement in TDOA estimation is a lower mean-squared error. Each method was tried 
using four different realizations for each of the three sampling times, for a total of twelve | 
different realizations per method used. However, we will only present one set from each 
sample time here, as each set follows a similar trend. The remainder of the plots can be 
found in Appendix [D]. To simplify explanation, we shall label the data sets as follows: 

a) HP-ADS generated signal with a sample interval of 92.59 nsec. 

b) HP-ADS generated signal with a sample interval of 185.185 nsec. 

c) HP-ADS generated signal with a sample interval of 370.370 nsec. 

d) Matlab generated signal. 
We will use the Matlab generated signal as our benchmark for comparison. For each 
method we will plot the mean-squared error against SNR’s in the range of — dB to 20 
dB. All Matlab codes used in these simulations are provided in Appendix (ЕЈ. For 
clarity in the following discussion, we define the terms average error and percentage 
improvement. Average error is calculated by adding all non-zero error values for a 
particular realization, and dividing that sum by the total number of the non-zero elements 
in that realization. Percentage improvement, as shown in the equation below, is 
calculated by subtracting the improved value from the original value and dividing the 
difference by the onginal value then multiplying by 100. 


original value — improved value . 
Да а SAO, percent improvement 


original value 
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І. 


Fourth Order Moment 


The results using this method are shown in Figure 5.4. There are two obvious 


trends that are noticed in this figure. We can easily see that the mean-squared error 


(a) 4th Order Moment HP-ADS; Ts = 92.59 nsec 
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(c) 4th Order Moment HP-ADS; Ts = 370.37 nsec 
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(b) 4th Order Moment HP-ADS; Ts = 185.185 nsec 
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(d) 4th Order Moment; Matlab 
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Figure 5.4: Results of Varying Sampling Interval for Fourth-Order Moment Technique 
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improves (decreases) as the sampling interval increases. This observation 15 presented 
numerically in Table 5.2. From this table, we can clearly see a decrease in error values as 
we move from left to right across each row. 

In the Matlab plot (Figure 5.4 (d)), or our benchmark plot, notice that the mean- 
squared error is zero at 6 dB and higher SNR for the Fourth-Order Moment denoised 


curve. Also notice that, as expected, the denoised curve is lower than the non-denoised 


Mean Squared Error 


case (a) case (b) case (c) case (d) 
Sample Time = Sample Time = | Sample Time = į Matlab Generated 
A поес 185.185 nsec 370.370 nsec 
No Fourth- No Fourth- No Fourth- No Fourth- 
Denoising | Order Denoising | Order Denoising | Order Denoising | Order 

Denoising Denoising Denoising Denoising 
26.30 11.04 5.80 
9.60 5.15 


3.82 2 lil 


2.39 





Table 5.2: Mean-Squared Error Values at Varying SNR for the Fourth Order Moment 
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curve at all points. This shows that the denoising is effective in producing an improved 
TDOA estimate. 

Using the HP-ADS data with a sampling interval of 92.59 nsec, case(a), the mean- 
square error at -3 dB SNR is 26.3 when fourth-order denoising is employed. This is a 
31.496 improvement from the non-denoised curve. The average percentage improvement 
over all SNR's between —3 and 20 1s 45.6%. 

Using HP-ADS data with a sampling interval of 185.185 nsec, we notice an 
average 52.296 increase with fourth-order moment denoising. The mean-squared error is 
5.68 at —3 dB, which is a 78.495 improvement over the —3 dB value of case (a). 

Using HP-ADS data with a sampling interval of 370.37 nsec, case (c), we obtain 
even better results. The mean-squared error at —3 dB is now 1.96, which is a 65.5% 
improvement over case (b). The average improvement realized by using fourth-order 
moment denoising instead of the non-denoised method in this case is 69.0%. We now 
see a definite trend that as our sampling interval increases, denoising increasingly 
improves the mean squared error values. 

In the benchmark case (d), there is a maximum mean-squared error of 0.25 at —3 
dB using fourth-order denoising. The average percentage improvement of fourth-order 
moment denoising over non-denoising is 72.5%. We also notice that the mean-squared 
error 1s zero for all SNR above 6 dB. 

The improvements we have discussed in this section can be related to the 
correlation function plots presented in Section 5.A. There, we saw that the correlation 


function plot improved with increasing sample time. Here, we state that a higher degree 


за 





mean square error 


mean square error 


of correlation of the signal components reduces the probability of error in the TDOA 
calculations. Thus, we achieve lower mean square error values with increasing degree of 


correlation of signal components. 


2. Modified Approximate Maximum Likelihood (MA ML) 


The results using this method are shown in Figure 5.5. We notice a similar trend 


as in the fourth order moment method. We end up with lower mean squared error values 


(a) MAML HP-ADS; Ts = 92.59 nsec (b) MAML HP-ADS; Ts = 185.185 nsec 
а ас лов са 10 T UN шыл Т н 
~ xcorr without denoising —~ xcorr without denoising 
Modified AML 9 4th order moment 
il 


1 7 
| 
| 






mean square error 
> ол 
— SE ES ee re a $1 05 


10 5 0 5 10 15 20 -10 5 0 5 10 15 20 
SNR (dB) SNR (dB) 
(c) MAML HP-ADS; Ts = 370.37 nsec (d) MAML, Matlab 
10 = ААА 


xcorr without denoising | 
| Modified AML | g| 





8 | 
7 | 7 
| 
| о 
6 = 6 
2 
5 a 5 
e 
c 
4 5 4 
E 
3 | 3| 
2 A 
1 T 
0 | 0c | 
-10 5 0 5 10 15 20 -10 -5 0 5 10 15 20 
SNR (dB) SNR (dB) 


Figure 5.5: Results of Varying Sample Times for MAML Technique 
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as the sampling interval increases. We also notice that these mean squared error values 
are lower than those obtained ın the fourth order moment method, as evidenced ın Table 
2 

In case (a), there is an average improvement in mean squared error of 89.6% by 


employing MAML denoising. Recall, from Section 5.B.1, the fourth order moment 


Mean Squared Error 
case (a) case (b) case (c) case (d) 
N {Sample Time = Sample Time = | Sample Time = | Matlab Generated 
92.59 nsec 185.185 nsec 370.370 nsec 


В [No МАМІ, [Мо МАМЕ | № МАМЕ № МАМЕ 
Denoising | Denoising | Denoising | Denoising | Denoising | Denoising | Denoising | Denoising 


11.04 las 0.63 0.09 
















0.01 






За 0.60 








0.27 







2:11 





0.13 


Table 5.2: Mean-Squared Error Values at Varying SNR for MAML Method 
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denoising method produced only an average 45.6% improvement for case (a). Thus, we 
have proved that the MAML denoising method is, in fact, superior to the fourth order 
moment. 

In cases (b) and (c), there ıs an average improvement of 89.5% and 88.8% 
respectively. These are in accord with the averages obtained in case (a). Thus, we have 
shown that through the use of the HP-ADS program, we can generate a simulated GSM 


signal and extract the necessary data to obtain a TDOA estimation. 


3. Time-Varying Approximate Maximum-Likelihood (TV-MAML) 

In Figure 4.1, we saw that the time-varyıng MAML (TV-MAML) method 
produced the best results at higher SNR. The data extracted from the HP-ADS program 
did not produce similar results as obtained in the fourth order moment technique and the 
MAML method. An unsuccessful run of the TV-MAML program is shown in Figure 5.6. 
We can see in this figure that the MAML plot outperforms the TV-MAML plot. 
Furthermore, the TV-MAML plot has a non-zero mean squared error even at very high 
SNR. This can be explained as follows. Recall from Section 4.C., the TV-MAML 
requires the MAML technique be applied twice to each segment of the signal ın order to 
reduce the probability of capturing non-signal components along with the actual signal. 
Recall also, all of our simulations used data sets that were 600 samples long. Therefore, 
the degree to which a particular data set 1s transitory depends upon its sample interval. 
Lower sampling intervals, as described in Figure 5.2, produce data sets that contain less 
information. Thus, the potential effectiveness of the TV-MAML is not fully realized if 


the signal does not contain a relatively large amount of information. 
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TV-AML Using HP-ADS Data Ts = 370.37 nsec 
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Figure 5.6: Results of Unsuccessful Run of TV-MAML Program Using Sample Interval 
of 370.37 nsec 

In order to evaluate the effectiveness of the TV-MAML method, we must use a 
higher sampling interval when collecting data from the signal. Thus, we found 740 
nanoseconds to be a suitable sample time. The signal ıs shown in Figure 5.7 along with 
the results from the TV-MAML program. Here, we see a definite improvement relative 
to Figure 5.6. The TV-MAML outperforms the MAML program at O dB (and lower 
SNR’s) by an average of 18.1%. This demonstrates that the TV-MAML works better in a 


higher noise environment than the MAML. 
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(a) HP-ADS signal | channel; Ts = 740 ns 
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Figure 5.7: (a) HP-ADS Produced GSM I-Channel Signal With Sample Interval of 740 


nsec. (b) Results of TV-MAML Using Input Signal of 740 nsec Sample Interval 
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VI. SUMMARY AND CONCLUSIONS 


A. SUMMARY 


The objectives of this thesis are to simulate GSM signals and to evaluate the mean 
squared-error for TDOA estimation. The Hewlett-Packard Advanced Design System is a 
very powerful communications development tool and proved invaluable toward our 
objective of GSM signal generation. Appendix [F] contains a concise user's guide to the 
HP-ADS program for performing the tasks conducted in this thesis. We were able to 
manipulate this system to provide signals with desired characteristics which would later 
be used in the Matlab environment to determine the effectiveness of three different 
denoising methods. The denoising methods were presented in earlier research. These are 
the fourth-order moment, modified approximate maximum-likelihood, and time-varying 
approximate maximum-likelihood methods. We used three signal sets generated by the 
HP-ADS program to test the performance of each of these methods. The results of these 
tests agreed with the earlier research. It showed that the MAML method is the best choice 
except when the signal contains a relatively large amount of information, in which case 
the TV-MAML tends to be better. 

In the process of testing the denoising techniques, we have demonstrated our 
ability to extract localization information from our generated GSM signals. This is 
manifested in the fact that we obtain relatively low mean squared error in our TDOA 
calculations. We also found that the error could be reduced further by increasing the 


interval over which samples are taken from the signal. 
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Although the Matlab programs provided the desired results, there were some 
inconsistencies to be noted. In the Matlab code for the TV-MAML, wavelet 
decomposition and re-composition 1s performed on the noise-infected input signal. At 
some point between the decomposition and re-composition, a dc offset is introduced as 
well as some slow frequency oscillation toward the end of the sequence. These issues 


will be addressed in subsequent work on this project. 
B. FUTURE WORK 


In addition to providing a solution to the oscillations and DC offset described in 
Section 6.A., it may be possible to modify the existing Matlab code to obtain even better 
mean squared error values. Follow on work should also re-examine the limitations of the 
time-varying approximate maximimum-likelihood method presented in Section 5.3. An 
extension to cases with different SNR’s is also needed to determine the performance in a 
fading environment. Also, the situation where there is a weak signal in one channel and a 
stronger signal in the other channel should be examined. This would provide insight into 
the question of what ranges of SNR are reliable. Finally, an investigation into the actual 
position location estimation using the solution of the hyperbolic equation presented in 


Chapter III should be conducted. 
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APPENDIX A. USA TODAY ARTICLE ON CELLULAR 911 LEGISLATION 


90% Washington 


10/26/99- Updated 05:39 PM ET 


Clinton signs nationwide 911 law 


WASHINGTON (AP) - President Clinton signed legislation Tuesday 
making 911 the official emergency number nationwide - for both 
regular and cellular phones. 

The measure also calls for development of technology that can track 
mobile callers. 


People with wireless phones now will be able to speed responses to 
highway accidents, crimes and natural disasters," Clinton said. 
"Getting rapid care tó someone who is suffering from a heart attack or 
Г involi i in a car crash can mean the difference between life and 

ea " қ 


While 911 is widely used as the emergency number for traditional 
phones, there are 20 different codes for wireless callers across the 
country. The changes are aimed at cutting response times for the 
emus who answer 98,000 emergency calls daily from cellular phone 
callers. 


"In my home state," said Sen. Conrad Burns, R-Mont., "three quarters 
of the deaths in rural areas are because the first responders couldn't 
get there in time." 


Health care professionals joined Burns at a Capitol Hill news 
conference to applaud the new law. 


"We have great emergency room personnel. We can do a lot for 
accident victims if we can find them and get them there," said 
Barbara Foley of the Emergency Nurses Association. "That's what 
this legislation helps us do." 


Another provision of the act directed the Federal Communications 
Commission to help states develop emergency systems, including 
technology that can automatically locate cellular callers who have 
dialed 911 or been involved in an accident. 


The FCC in September moved forward with plans to require that 
cellular 911 calls automatically provide a caller's location. Regulators 
want manufacturers to begin providing locator technology within two 
years. 
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Privacy advocates have raised concerns about potential abuse of the 
technology, which would take advantage of the Global Positioning 
System developed by the military. 


The law signed Tuesday called on regulators to establish "appropriate 
privacy protection for call location information,” including systems 
that provide automatic notification when a vehicle is involved in an 
accident. 


It said that calls could only be tracked in nonemergency situations if 
the subscriber had provided written approval. "The customer must 
erant Such authority expressly in advance of such use, disclosure or 
access." according to Senate documents detailing provisions of the 
legislation. 


An estimated 700 small and rural counties have no coordinated 
emergency service to call - even with traditional phones. The bill 
would encourage private 911 providers to move into those areas by 
granting the same liability protections to wireless operations that now 
are offered to wireline emergency service systems. 


Separately, the FCC took action earlier this year to increase the 
number of cellular calls to 911 that are successfully completed. The 
commission required that new analog cellular phones - not existing 
phones - be made with software that routes 911 calls to another 
carrier when a customer's own service cannot complete the call. 


Calls sometimes aren't completed because a caller is in an area where 
his or her carrier does not have an antenna, because networks are 
overloaded or because buildings or geography block signals. 


Digital phones, of which 18.8 million now are in use, were not 
covered by the new FCC rules adopted in May because such phones 
are more complex than their analog counterparts and there is no easy 
fix for the problem. 


#06 по: узт көч тарзе өө # e e e TTT 
290282554455 9504 68559 а ея ветвь дн ел ита aaa Вен оф вто о. баты о чо ety кат са 


e Go to Washington news 
e Go to News front page 





2 Сорупаћ: 1999 USA TODAY. a division of Gannett Co. inc. 
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APPENDIX B. HP-ADS GSM COMMUNICATION SYSTEM COMPONENTS 
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APPENDIX C. POWER SPECTRAL DENSITIES OF GSM TEST SIGNALS 


PSD COMPARISON FOR SAMPLE TIME = 92.592 nsec 
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PSD COMPARISON FOR SAMPLE TIME = 185.185 nsec 


Power Spectral Density Data Set 1; Ts = 185 nsec Power Spectral Density Data Set 2; Ts = 185 nsec 








Power Spectral Density Data Set 3; Ts 7 185 nsec Power Spectral Density Data Set 4; Ts = 185 nsec 
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PSD COMPARISON FOR SAMPLE TIME = 370.37 nsec 
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APPENDIX D. MATLAB SIMULATIONS FOR EACH DATA SET 


FOURTH ORDER MOMENT & MAML FOR SAMPLE TIME = 92.592 nsec 
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FOURTH ORDER MOMENT & MAML FOR SAMPLE TIME = 185.185 nsec 
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APPENDIX E. MATLAB CODES 


FOURTH-ORDER MOMENT DENOISING TECHNIQUE 


9. * * * * е Жк ke ck sk * * ke ke ke e e ke e he Se ec ke e ke ek ke ke ke ko e КЖ Кк Кк К A A ORO kO kO R K 
de c c Ж Ж Ж Ж Ж ж kk 
kkk kkk kkk kkk kkk kkk k kk kk k k kkk k k k k k k ККК Жо Жок К Ккк К Жк ЖКК К Ж k 


Ж Ж Ж Ж Ж Ж Ж Ж Ж Ж Ж Ж 


$ Denoise sta: Wavelet Denósing Based on The Fourth Order 


% SYNTAX: y=Denoise sta(xn, yn) 


RN DUI Zn Received signal from first receiver 
% уп = Received signal from second receiver 


% OUTPUT: у = Denoised signal Xn based on Үп statistics 


5 SUB FUNC; None 
Written by Spiros Mantis 


5 жжхххххххххххххххххххххххххххххххххххххххххххххххххххххххх 


oo 


+ + + + + + * + + * * + 
Qo ck ck * AA AAA AAA AAA AAA AAA AAA AAA AAA RAR 


ж ххх» * x 


function y=Denoise_sta(xn,yn); 


xyn=xcorr (xn,yn, 'biased'); 
[sigmas b]=max (xyn); 
rx=xcorr (xn, 'biased'); 
maXX=rX (length (xn)): 
ry=xcorr (yn, 'biased'); 
maxy=ry(length(yn)); 
sigmanl=maxx-sigmas; 
sigman2=maxy-sigmas; 


lamdax=3.1*sigmanl”2; 
lamday=3.1*sigman2%2; 


nx=floor (log2 (length (xn) )); 
ny=floor (log2 (length (yn))); 


Эй 


[cx lx]=wavedec(xn,nx,'db4'); 
[cy ly]=wavedec(yn,ny, 'db4'); 


ахе- Mr 

Тол ас Шарх 
d=detcoef (cx,1x,1); 
dl-length (d); 
A-(1/dl)*sum(d.^4); 


if A<lamdax 
dc=zeros(1,dl); 
else 
аста; 
епа 


dxc=[dc dxc]; 
end 


a=appcoef (cx, 1x, 'db4',nx); 
al=length (a); 
MEINES SUI A: 
if A«lamdax 
ac=zeros(l,al); 
else 
ac=a; 
end 
ерсе ао ерее: 


xd=waverec (dxc,1x,'db4'); 


y-xd; 
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© Ж * + зе e ke ck е k k k k k k k k k К ko ko kO k k ko SK ko Sk ko ko koko k S КК 
ххх Х 


Qo ck ck ck RAR RRA RARA + I * I T k k Tk k k k OK KOK OK kO kO K KOR KOK KOK KOK KOK e KOK KOK ke R KOK K 


dir dir dir di di di dir de +< + + + 


$ Modified tezó6: This is a test program for wavelet 


denosing 
$ based on the fourth order moment technique. 
5 GSM signal is used. 


ЕВЕ NIAX: моба тесш 26 
% INPUT: None 
$ OUTPUT: Mean square error versus SNR 


з» SUB FUNC: Denoise sta.m 


$ Written by Spiros Mantis 
Qu ck е е е е e eee k e e e e ee kO e e e e e e ee e eec e ee e e kO kO e ce e ke Se Ж RO R KO KO RO KOK O O RO O RO R 


жжжжжжжжжжжж 
Qs ck ke e e ee e he e k k e k Kk Tk k k K k ko T kO Ok e eee ce e e e eee he e eee e RO RO KR ORO R RO KOR OR KOR R 


Decke we e Je ce e eode e + + 


clear all 


oo 


gsm set; $Configuration variables created in memory. 
these are: 


oo 


5 Тр(е 3.692e-6) 

% ВТ(- 0.3) 

% OSR(= 4) 

° SEED(= 9515167065) 


TNTTEL = О 


оо 


oP оо 


s datafgen(INIT”"1I), 5 this creates a binary data 
(tx burst, 1,0]=gsm mod(Tb,OSR,BT,data, TRAINING); 


Ela; 

dataa O Serti 

s=transpose (s2); 
sl-length(s); 
pow-(1/sl)*sum(abs(s).^2); 


K=100 % number of realizations 
rand('seed',40); 
A Ian а k); 
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delay=floor (f); 
delay(1:K/2)=-delay(1:K/2); 5 delay is between -150 to +150 


n= [1207782177 1615 14 ИСЗ 0 366 
БМВ щи): 
гора пе шешо сие.) 


oran=SNR(k) 
for i=1:K 


x=[zeros (1, 200 11582805 (1,224)1; 
y=|zeros(1,2001G b y i)) s zeros(@ 224 delay (1))1% 


randn('!'srater Бу); 
noil opeal-sghEbRpew/ (2*oran))*rameuw ІШЕДІ); 


randn ('state',3*(i+j)); 
noil imag=sqrt (pow/ (2*oran)) *randn (1,1024); 


тапап (“бәс (its) ); 
noi2_real=sqrt (pow/ (2*oran) ) *randn(1,1024); 


randn('state',5*(1i+3)); 
noi2 imag=sqrt (pow/ (2*oran))*randn (1,1024); 


nomas ea noil imag; 
nouwecnor2resltj*noi2 imag; 


xn=x+no11; % 


x 
yn=y+no12; % y + noise 


% IDOA caleulation with xcorr( without denoising) 


о 


xy=Xcorr(xn,yn); + correlation о: = апас таа ее 
[al b1]=max (real (xy)); 


ег1 (1) =Ае]ау (1) - (61-1024); 


2 Fourth order moment denoise 
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X real=Denoise sta(real(xn),real(yn)); 
X imag=Denoise sta(imag(xn),imag(yn) ); 
Yrreal-Denarsenst ac (vn),reallzn)); 
Y imag-Denoise sta(imag(yn),imag(xn)); 


X-X real*j*X imag; 
Y-Y real*j*Y imag; 
SACO ОС 


$ Correlation of X and Y (denoised) 
[a2 b2]=max (real (XY)) 


е 
7 


erl0(i)=qelay(1)À- (2 ` 1024); 


end 


errorl0a370set4(k 


)=(]Z]enoth l er10))*sum(er10.22); 
errorlsta370set4(k)= 


tenia (ex) )*sum(erl.*2); 


H10a370set4(k,:)=er10; 
Hsta370set4(k,:)=er1; 


end 


Елаике (6) 

ВО 1898 17 lG 15 14 13 12 9 6 3.0 -3 -6]; 

plot (k,errorlsta3/0Oset4 (1:14), 'o',k,errorl0ßa370set4 (1:14), ' 
x',k,errorlsta370set4 (1:14) ,k,errorl0a370set4(1:14)) 
legend('xcorr without denoising','sta') 

title('4TH ORDER method with cross corellation after 
semoise; T00 realizations of data 370 set4') 

ylabel ('MSE') 

xlabel('SNR') 

figure (7) 

plot (1:2047,xvy) 

title('Correlation Function of x and y signals w/noise') 
figure(8) 

lot (1: 2047 XY) 

title('Correlation Function of x and y signals without 
noise') 
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Save 
Save 


save 
Save 


errorj0a370set4; 
errorlsta370set4; 


Hsta370set4; 
H10a370set4; 
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MODIFIED APPROXIMATE MAXIMUM LIKELIHOOD DENOISING 
TECHNIQUE 


EP ck ck Ж З З e * ke e * e * k he ke se k е k k k e e he k k kok ko ko koko o k k k o ecce ko koko k kG ko Ó kG К k Sk Ж Ж 
Ж ж + + Ж Ж Ж Ж Ж Ж Ж Ж 

% * * * ХХ e ke se k k k k k k k k kO k e ke ce eee e e ke e e e ck ek e k k ko ke ke ke A ke Хк Ж k SK 
Ж Ж Ж Ж Ж Ж Ж Ж Ж Ж ж ж 

Denoise: Approximate maximum-likelihood delay estimation 


via orthogonal wavelet transform. In this function, 
we assumed that noise is Gaussian noise and it has 


flat freq response. We modify each detail function 
by multipliying modified AML coefficient based on 
the signal and noise powers. 


SYNTAX: y-amll (xn, yn) 


INPUT: xn = Received signal from first receiver 
yn Received signal from second receiver 


OUTPUT: 
y = Denoised signal 


SUB” EUNC: None 
Written by Spiros Mantis 


хх ЕХ ХКК ХКК A A A AI A A A A AI Ж 


оо до до до до до до до до оо до до до о0о до (Ј 00 оо оо 


o 


5 
Ak A A AA RA Ax A A ЖК 
з. ck ck ck A + ke hee ХХ ЕХ КК kok k k K ko k kok ko ko ko Sk ko ko ko kk ko ko ko 


Ak + + + * k + * A AR 


function y=denoise (xn, yn); 


xyn=xcorr (xn, yn, 'biased'); 
[sigmas Юр] =тах (хуп); 
rx=xcorr (xn, 'biased') 
maxx=rx (length (хп)); 
ry=xcorr (yn, 'biased') 
maxXy=Yy (length (yn) ); 
sigmanl=maxx-sigmas; 
sigman2=maxy-sigmas; 


“е 


“е 
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nx=floor (log2 (length (xn))); 
ny=floor (log2 (length (yn))); 


[cx 1x]=wavedec (xn,nx, 'db4'); 
[cy ly]=wavedec(yn,ny,'db4'); 


dzem 

foul: nx 
d=-deseeetiex,1x,i); 
dl-length (d); 
sigmad- (1/dl)*sum(d.^2); 
Sigmasd-sigmad-sigmanl; 


sigmaa (i)=(2* (nx-1)) *sigmasd; 
if sigmasd<=0 

wd=0; 
else 


wd=sigmasd/ (sigmanl*sigman2+sigmasd* (sigmanl+sigman2)); 
end 
dc=wd*d; 
(C = [de dxc]; 

end 


a=appcoef (cx,1x, 'db4',nx); 
al=length(a); 
aux=sum(sigmaa) 


sigmasa-((2^nx))*sigmas-aux; 
if sigmasa<=0 

wa=0; 
else 


wa=sigmasa/ (sigmanl*sigman2+tsigmasa* (sigmanl+tsigman2)); 
end 
ac=wa*a; 
dxc=[ac dxc]; 


xd=waverec (dxc, 1x, 'db4'); 
y=xd; 
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Ss Sk ck k e e Se ke e he he e ee e ke ee he ke e he e e е TH TR TH TFT TI TI TITTEN N 


ххжххххххххх 


Q k k k k Ж Ж 


хжхххххххххх 


со 


5 Мосіштес рез о Y P Es St program for moq ri JR N. 
technique. In this 
program, we used GSM signal. 


до 90 


до 


о МА МОС Ене ez, 


oo 


ою 


INPUT: None 


ою 


oo 


QUITPUTS Mean square error versus SNR 


ою 


ою 


SUB ОМ ВЕБ. Ш 
Written by Spiros Mantis 


схххххххххххххххххххххххххххххххххххххххххххххххххххххххххх 


ою 


WX ok Wo kk otc ke ox xà x 
C ke ke ck coke ke ke ke ke e ke e ke e e eee e ke e ke ХКК ХХХ e ke К e e e e ke ke ke à à e e à à A + Ж Ж 


хххххххххххх 


clear all 


gsm set; $ Configuration variables created in memory. 
$these are: 


ою 


Tb(= 3.692e-6) 
вее) 

OSR(= 4) 

SEED(= 931316785) 
парт ce 


до до 


до 


ою 


data-data gen(INIT L); $ this creates a binary data 
[tx burst,I,O]-gsm mod(Tb,OSR,BT,data,TRAINING); 


SOS 
51=ЈепаЕћ (5); 
ром= (1/51) *ѕит(арѕ (5). ^2); 


K=100 $ number of realizations 
rand ('seed', 40); 
f=150*rand(1, Ky; 
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delay=floor (f); 
delay (1:K/2)=-delay(1:K/2); 5 delay is between -150 to +150 


п= | СЕС пенса 14-13 127905 3 0-50-61; 
SNRS 1052 (r. 0 
for k=1:length (SNR) 


oran=SNR(k) 
БӘ ІЗ 


x=[zeros (1,200) s zeros(1,224)]; 
y=[zeros (1,200+delay(i)) s zeros(1,224-delay(i))]; 


анал! state 2" (1, 
noil_real=sqrt (pow/ (2*oran) ) *randn (1,1024); 


randn('state',3*(1+3)); 
noil_imag=sqrt (pow/ (2*oran) ) *randn (1,1024); 


randn('state',4*(i*j)); 
noi2_real=sqrt (pow/ (2*oran) ) *randn (1,1024); 


и-шеши искате! БН (11) 12 
noi2 imag=sqrt (pow/ (2*oran))*randn (1,1024); 


noil=noil_real+j*noil imag; 
noi2-noi2 real*j*noi2 imag; 


xn-x4noil; $ x + noise 
yn=ytnoi2; % y + noise 


$ TDOA calculation with xcorr( without denoising) 

Xy-cxcorrdixm vH S correlation orf Tande With хобі 
$(with the presence of noise) 

[al bl]=max (real (xy)); 

erl(i)=delay(i)-(b1-1024); 


SAML 


x real=denoise(real(xn),real(yn)); % Denoising of the 
received signals 
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Ишиса пепел ве (ези, веат); 
x imag=denoise (imag (xn) ,imag(yn)); 
y imag=denoise (imag (yn), imag (xn)); 


X=x real+3*x imag; 
Y-y real*j*y imag; 


KI-ZeorriX 6) > 


$ Correlation of X and Y (denoised) 
[a2 b2]=max (real (XY) ) 


er8 (i)-delay(i)-(b2-1024); 
end 


error8a(k)-(1/length(er8))*sum(er8.^2); 
errorla(k)-(1/length(er1l))*sum(er1.^2); 


H8a(k,:)=er8; 
Hla(k,:)=erl; 


end 


figure(6) 

Е- |20 18 1 ТЕ ІС ТА 13 12 9 6 3 0 3261; 

puro E, cM Га MGUUTEUerroresswd1:14),;'x'-k,errorla(l:l14 
ПАК етвоква (1: 14)) 

legend('xcorr without denoising','aml') 

title('Denoise AML method with cross corellation after 
denoise; 100 realizations') 

ylabel ('MSE') 

xXlabel ('SNR') 

figure(7) 

Ісі (1:2047, ху) 

title('Correlation Function of x and y signals w/noise') 
figure (8) 

роі 294 7, ХУ) 

title('Correlation Function of x and y signals without 
noise') 


save errorßa; 
save errorla; 


save Hla; 
save Н8а; 
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TIME-VARYING MODIFIED APPROXIMATE MAXIMUM LIKELIHOOD 
DENOISING TECHNIQUE 


L k kk k k хк хк хх Ñ Tc TH Ñ k жск ху Ку к ККК к хх КККК ko ko ko kO ko ko Sk ko k ko K k a Sa 

ЖЖ КОЖ X W W < X X x X 

Я ke у» +++» +++ ХХХ ke ko ko Sk Sk Sk k Sk Sk k ko Sk ko k К К k k e k e К K SK SK 

хх Хе 

aml2: We modified the MAML technique by dividing each 
detail function into two segments. Different 
coefficients for each segment are computed. 


до до с о 


ою 


SYNTAX: y=aml2 (xn, yn, delay) 


ою 


$ INPUT: xn = Received signal from first receiver 
5 уп = Received signal from second receiver 
5 delay = True TDOA between xn and уп 


ою 


OUTPUT: y = Error between true TDOA and estimated TDOA 


ою 


ою 


SUB FUNC: None 
Written by Unal Aktas 


ххх ХКК k K e kk ke e ke e e e ee e kok e e e e e e e e e e e e e Ж 


ою 


жжжжжжжжжжжж 
ЖЖЖЖ АККА АССА СААС АССА ССА АССА АССА А АС С АС А С З С С К 


dee ee se ee Ke ee < 


function y=aml2 (xn, yn,delay) ; 


xyn=xcorr (xn, yn, ‘biased ); 
[sigmas b]=max(xyn) ; 
rx=xcorr (xn, 'biased'); 
maxx-rx(length(xn)); 
ry=xcorr(yn, 'biased'); 
maxy=ry (length(yn) ); 
Sigmanl=maxx-sigmas; 
sigman2=maxy-sigmas; 


nx=f£floor (log2 (length (xn))); 510 
ny=floor (log2 (length (yn))); 
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[cx 1x]=wavedec (xn,nx, 'db4');%3wavelet decomposition of xn 
Sat level 10 
(су ly]=wavedec (yn,ny, 'db4'); 


Єхс=|]; 

Па о а а — |: 1 
dedetcoef(cx,lx,i); $detail coefficients at level i 
dl=length(d); 57 
№51=128/2^ (1-1); % the length of the subblock 

Tf NSI«-1 912-09. 1<7> 0 


Ns-dl; $7 
else 

Ns=Ns1; % .25<М51>128 
епа 


D=ceil(dl/Ns); 
if dli«Ns*D 
dm=[d zeros(1,D*Ns-dl)]; 
end 
for k=1:D 
p=(k-1)*Ns+1:k*Ns; 
sigmad=(1/Ns) *sum(dm(p) .*2); 
sigmasd=sigmad-sigmanl; 
if sigmasd<=0 
wd=0; 
else 


wd=sigmasd/ (sigmanl*sigman2+sigmasd* (sigmanl+sigman2)); 


end 

qe (p) маецца 5 
епа 
ахе=гае( 1 ге аа: 


епа 


a=appcoef (cx,1x, 'db4',nx); 
al=length (a); 
sigmaa-(1/al)*sum(a.^2); 
Sigmasa-sigmaa-sigmanl; 
if sigmasa<=0 

wa=0; 
else 

wa=sigmasa/ (sigmanl*sigman2+tsigmasa* (sigmanl+sigman2)); 
end 
ac=wa*a; 
ахс= [ас ахс]; 
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xd=waverec (axe, lx, dol 1; 


dyc=[]; 

for i=l my 
dy=detcoe (Cy, ly, i); 
dyl=length (dy); 


Ns1=128/2* (1-1) ; $ the length of the subblock 
1f Nsi<=1 
Ns=dyl; 
else 
Ns=NS1; 
end 


D-ceil(dyl/Ns); 
неу 1 <Ns>D 
dmy=[dy zeros (1,D*Ns)]; 
end 
for k=1:D 
p=(k-1)*Ns+1:k*Ns; 
sigmady-(1/Ns)*sum(dmy(p).^2); 
sigmasdy=sigmady-sigmanl; 
if sigmasdy<=0; 
wdy=0; 
else 


wdy=sigmasdy/ (sigmanl*sigman2+tsigmasdy* (sigmanl+sigman2)); 
end 
dcy (p) =wdy*dmy (p) ; 
end 


8  [dcy (1:dv1) dye]; 
end 


ay-appcoef (cy,ly,'db4',ny); 
ayl=length (ay); 
sigmaay-(1/ayl)*sum(ay.^2); 
sigmasay-sigmaay-sigmanl; 
if sigmasay<=0 

way=0; 
else 


way=sigmasay/ (sigmanl*sigman2+sigmasay* (sigmanl+sigman2)); 


end 
acy-way*ay; 
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dyc=[acy dyc]; 
yd=waverec(dyc,ly, 'db4'); 


$5 [ма1 ма2 маз ма4 ма5 маб ма7 мав ма9 мс9) 
аа Corr (xd ya 


Га5 b5]=max (rxyd) ; 


er5=delay-(b5-1024) ; 
угеко; 


A 


QS ck kk * * * * * T * * * * * * e e e ke e e e he e ee he he e e ehe TI TH ce ж ж жж ж A A A A A К ke ke A A Ж Ж Ж Ж Ж 
ххх жя 

ЖКЖ ЖЖЖЖ ЖЖК ЖК К ЖКК К К Ж Ж К К К К К К К А А А К с К А К с З с А дс С К А д С К С ЗК 
хххххххххххх 

[e] е Ц . . а 

6 tez5t : This is a test program for time varying MAML 


technique. 
In this program, GSM signal is used 


оо 


оо 


SYNTAX: tezb5t 


до 90 


до 


INPUT: None 


oo 


оо 


ние can squarezerrorzversuszshnk 


o9 


оо 


эра ые ат 2. 
Written by Unal Aktas 


ХК ХК КАЖЕ Ж +Х 


оо 


хххххххххххх 
хх eee e e ee e ee ee ee eee e e e ecc eee e ee ecce e e e ee e e e ko Ук Жу à Ax 


хххххххххххх 


clear all 


gsm сев; ЕЕ Lion varr Dic с леска сесии пеели ша 
шесе аге: 


ою 


ТЬ (= 3.692е-6) 
Bo 

OSR(= 4) 

SEED(= 931316785) 
IN O 


со до до 


ою 


Са а“ ава уеп (ГЫ И о Еле а ен се = а о а ан Ma 
[tx burst,1I,Q]=gsm_mod (Tb, OSR, BT, data, TRAINING) ; 


$s-Itj*Q; 

data TONS CI) 
s-transpose(s2) 

sl=length (s2); 

pow- (1/sl)*sum(abs(s).^2); 


K=100 % number of realizations 
rand('seed',40); 
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ре Бан K; 

delay=floor (f); | 
delay(1:K/2)=-delay (1:K/2); 5 delay is between -150 to 
+150 


ASIA о G №5 14 13 ЗО |, 
ӘМЕСІС УТО): 
for k=-Palengeh‘SNk) 


oran=SNR(k) 
hor К 


х-риснез (1259) = ?гегов (1,224)); 
v [2eros (1, 259+delay гето Сге (пр то 


randn('state',2*(i+j)); 
noil_real=sqrt (pow/ (2*oran) ) *randn (1,1024); 


ranenkesstate’,3*(i+4j)); 
noil imag=sqrt (pow/ (2*oran) )*randn (1,1024); 


randn ('state',4* (i+j)); 
ПЕ ЕЕ Боарон ози randn(1, 1024); 


тапай вгавс , 59 (141) ре 
noi2 imag-sqrt (ром/ (2*огап)) *гапап (1,1024); 


noil=noil_real+j*noil imag; 
nol2 msn roll плот 2 пад; 


xn=x+noil; % x 
yn=y+no12; 5 y + noise 
ADD Ea ELO with со without denoising) 


УЕ CO ny) o Correlation of x and y with xcorr 
[al bl]=max(real(xy)); 


erl (1) =delay(1)-(b1-1024); 


SAML 
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el=aml2 (real (xn),real(yn),delay(i)); 
e2=aml12 (imag (xn) ,imag (yn) ,‚delay(i)); 


ela=aml1 (real (xn), real (yn) delay(i)); 
eZ2a=aml1 (imag (xn), imag (yn) delay(i)); 


er8(i)=(el+e2)/2; 
ег8а (і) = (е1а+е2а) /2; 
епа 


error8t(k)-(1/length(er8))*sum(er8.^2); 
errorla(k)-2(1/length(erl))*sum(er1.^2); 


error8a(k)-(1/length(er8a))*sum(er8a.^2); 


H8t(k,:)-er8; 
Hla(k,:)=erl; 


end 
$load error8a 


figure(6) 
О ІТ Ир е еды 
Bor K,errorlal1l214)72 075, K, errorse | 4s x error Ban] Ta 
INCOME 
k, errorla({1:14) Kk error8t(1:14),;,k errorga(1:14)) 


legend('xcorr without denoising','time varying aml','aml') 
title('TVAML Method ADS input data; Ts = 740nsec Set#1; 100 
realizations') 

figure(7) 

plot(1:541,real(s2)) 

title('ADS signal (s) real component; 740ns') 


save error8t; 
Save errorla; 


save Hla; 
save H8t; 
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APPENDIX F. GUIDE FOR USING GSM SIMULATION IN HEWLETT- 


PACKARD ADVANCED DESIGN SYSTEM 


STEP 1: Start Program 


If there is no icon on startup window: 
Start => Programs => HP-ADS => ADS 


STEP 2: Open GSM Example Project 


From ADS (main) window: 

File => Example Project => Com_Sys => gsm_prj => ok => wait a few seconds 

for windows to pop up. (two windows come up: GSM_SYS Schematic:2 and MODEM 
Schematic:]. The Modem Schematic is for a more detailed simulation. I have been primarily 
working with the GSM_SYS schematic.) DO NOT MAKE ANY CHANGES TO EITHER OF 
THESE FILES!!!! Any changes will be saved in the system and they will remain there for future 
users of the program. Save to your H: drive or other drive of your choice and make changes there. 


STEP 3: Modifying the system 
The options here are unbounded. You can experiment with all the different components in the 
Component Library. But for now, we just want to know how to get data out of this system. 


This is done by placing ‘data sinks’ in the desired locations. 


Go to the common component box at the top of the schematic window and click the 
arrow. This will show all the different components available. => 


For data extraction, you want to use ‘Sinks.’ From there you can find the appropriate 
‘sink’ for your data. (Timed and Spect. Analyzer are the most commonly used)=> 


Click once on the sink of your choice => 


Move the cursor into the schematic area and place the sink in the desired position by left- 
clicking when you have it in a clear spot. => 


End the command by using the right-click button on the mouse then click ‘end 
command’ or just hit the escape key. => 


Use wire to connect it to the system. (there is a wire button on the third row of toolbar, or 
you can go to the component button at the top) => 


Notice the name of the sink. (this is listed under the device type) This is important, 
because you will need to recall this name when you want to plot or export the data. You 
can change the name by clicking once on it and then typing the new name right there in 
the same spot. 
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STEP 4: Simulations 


After you have your desired system configuration, are ready to simulate the system: 
Simulate(on top toolbar) => simulate => wait for “simulation finished” in summary window. The 
simulation is now complete and you can extract your data. 


STEP 5: Displaying Data 


Window(on top toolbar) => New Data Display => wait for data display window 

To see an x-y plot: 
Click once on the “plot” icon => move cursor into white area of display(a box appears), 
position the box and left-click. => choose what data you want to plot. => ok => your plot 
appears in the data display window. 


To see data: 
Click once on the “123456” icon => follow the same procedure as above. 


STEP 6: Extracting and Exporting Data 
I] haven't figured out how to export the plots but here is how to extract the data: 
In the schematic window: click on "start the instrument server” button(third button third row of 
toolbar) => Click the following: 
Read/Write: Write 
Write to: File 
File Format: Citifile => 
Use the browse button to choose a directory and file name. => 
Choose desired data from Datasets box. => 
Click ‘write to file’ button. => your data is now saved in the specified directory. 
You will have to use wordpad or notepad to view this data. Then you will have to sort through all 
the data given to get to the specific data. This may be annoying, but rest assured the data you want 
is in there somewhere!! 
STEP 7: Importing your ADS data into Matlab 
After you have located your data in wordpad you can cut and paste from there into Matlab. 
If you want to know how many data points you have it is best to use the Matlab Editor/Debugger 
(unless you have a few hours to spend counting). 
- Make the data into a vector in Matlab and use the size command to determine the 


length. 
- Using this method you can cut it down to the right size. 
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